An order parameter equation for the dynamic yield 
stress in dense colloidal suspensions 



Michio Otsuki and Shin-ichi Sasa 

Department of Pure and Applied Sciences, University of Tokyo, Komaba, Tokyo 
153-8902, Japan 

E-mail: otsukiOj iro . c.u-tokyo . ac. jp, sasaOjiro.c.u-tokyo. ac.jp 

Abstract. We study the dynamic yield stress in dense colloidal suspensions by 
analyzing the time evolution of the pair distribution function for colloidal particles 
interacting through a Lennard- Jones potential. We find that the equilibrium pair 
distribution function is unstable with respect to a certain anisotropic perturbation in 
the regime of low temperature and high density. By applying a bifurcation analysis to 
a system near the critical state at which the stability changes, we derive an amplitude 
equation for the critical mode. This equation is analogous to order parameter equations 
used to describe phase transitions. It is found that this amplitude equation describes 
the appearance of the dynamic yield stress, and it gives a value of 2/3 for the shear 
thinning exponent. This value is related to 6 in the Ising model. 
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1. Introduction 

Soft glassy materials, such as dense colloidal suspensions and super-cooled liquids, 
exhibit diverse rheo logical phenomena ^1 12] ■ Typical examples are the decrease of the 
viscosity with increasing the shear rate, called shear thinning, and the increase of the 
viscosity with the shear rate, called shear thickening. These phenomena appear when 
the shear rate is beyond the Newtonian regime, that is, the regime in which the shear 
stress axy depends linearly on the shear rate 7, and hence the viscosity, rj = cr^y/'j, is 
independent of the shear rate. Interestingly, the Newtonian regime becomes narrower 
as the temperature decreases and the density increases, and eventually it disappears. In 
the case that there exists no Newtonian regime, the stress is finite in the limit 7 \ 0. 
This stress is called the dynamic yield stress. We seek to understand these rheological 
properties systematically on the basis of a microscopic description of the system. 

Recently, the nonlinear rheology of soft glassy materials has been studied 
extensively through use of molecular dynamics simulations IH E] , analysis of random 
spin models [HllZl, and mode coupling theory applied to systems under shear [Hj lU} ITU]. 
In particular, the mode coupling theory have predicted that the dynamic yield stress 
appears discontinuously at the glass transition point [H] and the results have been 
compared with a numerical simulation [TT| and an experiment ^21- On the other hand, 
the power law r] ~ 7"^/^ has been observed for a shear thinning fluids IE]- 

Among theoretical approaches, one standard method for microscopic study of 
nonlinear rheology is based on the mode coupling theory. In this theory, the singular 
behavior of the viscosity is described by the anomalous part of the time correlation 
function, using a generalized Green-Kubo formula JH] • As another approach describing 
the singular behavior, one could focus on the shear stress. Here, let us recall that the 
shear stress a^y is defined as the average of the y component of the inter-particle force 
acting on a plane transverse to the x direction. When there exist only two-body forces 
among the particles, this average can be expressed in terms of the pair distribution 
function. Therefore, in this case, the rheological properties mentioned above can be 
accounted for through an analysis of the pair distribution function. 

Employing an approach of the latter type, in the present Letter, we first investigate 
the linear stability of the equilibrium pair distribution function. We find that this 
function is unstable in the regime of low temperature and high density. Next, applying 
a bifurcation analysis to a system near a critical state at which the stability changes, we 
derive an order parameter equation that describes the appearance of a dynamic yield 
stress in a simple manner. 

2. Model 

We consider a system consisting of spherical colloidal particles suspended in a solvent 
under shear flow described by v{r) = (7?/, 0,0). We denote the volume of the system 
by V and the temperature by T. Let F = (ri, ■ ■ ■ , ttv) represent the particle positions. 
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The time dependent distribution function for F, \l/7v(r,t), satisfies the Smoluchowski 
equation pP 

= Ev. • (|v. - if.(r) - *.(r,i). (1) 



Here, the Boltzmann constant is set to unity, R is the friction coefficient, and -F'i(r) is 
the interaction force defined by Fi(T) = — Vj Iii this work, we employ 

the Lennard- Jones potential, with the explicit form V{r) = 4e (r/a) — (r/a) ^ . In 
the treatment below, a, e and R are set to unity, and all quantities are converted to 
dimensionless forms. 

From the N particle distribution function \l/jv(r,t), the pair distribution function 
g{r,t) is defined as 

gir, t) = V'J Sr.d^n ■ ■ ■ d^r^^^iT, t), (2) 

where r = ri — r2. Using this function, we can express the time-dependent shear stress 
as dH 

..,w4/dV«(M)^^^, (3) 

where p = N/V is the average number density of the colloidal particles, r = {x,y,z), 
and r = \r\. Here, we have ignored the hydrodynamic contribution to the shear stress. 
It is useful to expand g{r,t) in spherical harmonics, and we write 

g{r, t) = fj{r, t)lmY2,2{9, 0) + /«(r, t)ReY2,2{9, 0) 

+ J2 Gi,mir,m,Ue,(l)), (4) 

l>0;\m\<i,{l,m}^{2,±2} 

employing the spherical coordinate system {r,9,(j)). Hence, substituting equation (jH) 
into equation (jH}, we obtain 

a.yit) = ^p^fdrr^^f,ir,t). (5) 

This expression clearly indicates that if there is singular behavior of the shear stress 
ar^y, that of the pair distribution function can be observed. In particular, the existence 
of a dynamic yield stress implies that the pair distribution function in the limit 7 \ 
differs from the equilibrium one. This observation naturally leads us to conjecture that 
the dynamic yield stress is related to the instability of the equilibrium pair distribution 
function, gcq{r). For this reason, we carry out a linear stability analysis of geqij")- 

In order to determine the linear stability of 5'eq('"), we need to study the time 
evolution of the pair distribution function. This is described by the BBGKY hierarchy, 
and hence it is determined by the three-particle distribution function g^{r,r',t). 
In order to obtain a self-contained description, we truncate the BBGKY hierarchy 
by employing the Kirkwood superposition approximation, assuming the relation 
(73(r,r',t) = g{r,t)g{r' ,t)g{r — r',t) This approximation has been used in the 

calculation of the pair distribution function for both an equilibrium system ^E] and 
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a non-equilibrium system under shear P^l- With this approximation, the evolution 
equation for the pair distribution function is derived as 

dg{r,t) 



at 

with 



-V-J(r,t), (6) 



J(r,t)= -2T\/g{r,t)-2VVir)gir,t) 

- 2p i^J dV'Vy(r')(7(r, t)g{r', t)g{r - r\ t) 

+ iy^g{r,t). (7) 
3. Linear stability analysis 

First, the equilibrium pair distribution function geqix) for this model is obtained as the 
isotropic solution of the equation J = in equation ((Tj) with 7 = 0. This equation 
is called the Born-Green equation jl5j, which has been solved numerically ^Hj- Then, 
writing 

g{r,t) = g,^{r){l + h{r,t)), (8) 
we substitute equation (jSJ into equations (0) and ^ and obtain the form 

^eq(r)^/i(r, t) = C{h{r, t)) + M{h{r, t)) + ^g{h{r, t)), (9) 

where £ is a linear operator and M contains only second and third order polynomials 
in h{r, t). 

Because the equilibrium pair distribution function geq{r) is non-negative, it is 
linearly unstable if and only if the linear operator C has a positive eigenvalue. Therefore, 
in order to determine the linear stability of geq{r), we could numerically compute the 
eigenvalues of C. However, due to the three-dimensional spatial dependence of h{r,t), 
the computation of these eigenvalues is not simple. Therefore, to simplify the problem, 
we assume perturbations h{r,t) of the form h{r,t) = '?/'(r, t)ImF2,2(^, 0)- Such an 
assumption is reasonable, because the equilibrium pair distribution function is expected 
to be unstable with respect to perturbations of the form Imy2,2(6',0) in the regime 
characterized by low temperature and high density. (See equations (jU and (0)). Then, 
using the explicit form of C, we can rewrite C{h{r,t)) in equation as 

C{h{r, t)) = lmY2,2{e, <j))M{^{r, t)). (10) 

To solve the linear stability problem, we seek the eigenvalues of Ai. 

The eigenvalues of the operator Ai are computed numerically in the following way. 
We first restrict the spatial domain of ip{r,t) to the interval [0,1], with the boundary 
conditions dipir, t)/dr = a.t r = and r = I. We next approximate the linear operator 
as a matrix M by using a difference method with a spatial mesh size Sx ^HI- We 
calculated the eigenvalues of M for many values of {p,T). 
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In figure ^ we display the stability diagram obtained numerically. It is seen that 
the pair distribution function Qeqir) is unstable in the low temperature, high density 
regime. Note that the form of the boundary between the two regimes is qualitatively 
similar to the curve representing the glass transition in the (p, T) plane calculated using 
mode coupling theory jTHj, though at each density, the value of the temperature on 
the boundary in the present case is approximately twice that in the glass transition 
case. This overestimate of the temperature might be caused by the inaccuracy of the 
equilibrium pair distribution function geq{r) calculated with the Kirkwood superposition 
approximation. Indeed, it has been known that the maximum of gcq{i^) obtained by this 
approximation is less pronounced and is shifted towards smaller interparticle distances 
than that obtained by Monte Carlo simulations. 



3 

T 



2 M 

1 1.04 p 1.08 1.12 

Figure 1. Stability diagram in the {p,T) plane. The squares indicate states where 
the maximum eigenvalue of M is positive, and hence where 5cq('") is unstable. This 
result was obtained for the numerical values I = 7.0 and Sx — 7.0/512. 



4. Non-linear analysis 

Next, we focus on systems near the critical state at which the stability changes. If we 
fix the density, then we find that there exists a critical temperature Ts below which 
the equilibrium pair distribution function geq{r) is unstable. Let be the critical 

eigenfunction of the operator Ai at T = Tg. Then, we write a perturbation h{r,t) in 
the form 

h{r, t) = A(t)^,(r)Imr2,2(^, 0) + s(r, A{t)). (11) 

Here, s(r, A{t)) represents the contribution to /i(r, t) that is not from the critical mode, 
and we have assumed that its time dependence is restricted to that of the amphtude 
A. This is a reasonable assumption because the amplitudes of non-critical modes decay 
quickly to the values determined by A PU] . 

Then, using a bifurcation analysis (See reference |25 ^is a review), from equations 
and (fTT|) . we perturbatively obtain the evolution equation for A{t) in the form 

^Mt) = (Ts - T)aA{t) + G{A{t)) + iH{A{t)), (12) 
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where a is a positive constant, and the functions G{A) and H{A) can be calculated 
from equation [22] • In this Letter, we do not calculate G{A) and H{A). Instead, 
considering the general forms of these functions, we qualitatively describe the singular 
behavior of the shear stress. 

First, because the system possesses symmetry with respect to the simultaneous 
transformation x — > —x and 7 —>■ —7, G{A) must be odd function and H[A) an even 
function. We therefore expand as 

G{A) = hA^ + hA'' + ■ ■ ■ , (13) 

H{A) = Co + C2A^ + ■■■. (14) 

Next, let A^..{'~f) be a stable stationary solution of equation p2j). Then, from 
equations © and (fTT]) . the shear stress for systems near the critical state T = Tg 
and 7 = can be approximately expressed as 

[2% roo dUir) 
a^yC^A,{^)^—p^j^ drr^—^g,^{r)tlj,{r), (15) 

because s{r,A{t)) in equation (|TT| contains no terms linear in A as we have confirmed. 
Thus, the value of ^4^(7) determines the behavior of the shear stress near the critical 
state. In this sense, the amplitude of the critical mode can be regarded as an order 
parameter describing the appearance of the dynamic yield stress. 

Now, we investigate equation (fT^ . For 7 = 0, this equation possesses the trivial 
stationary solution A^ = 0. This solution is stable only when T > Tg. The non-trivial 
stationary solutions depend on the form of G{A). First, let us consider the case that 
63 < in equation (fT^ . (We conjecture that 63 < for the model we study.) Then, for 
T ~ Ts and 7 ~ 0, A* satisfies (Tg — T)aA* + b^Al ~ — 7C0 under the assumption that 
A^ and 7 are expressed as power-law functions of (T — Tg). Thus, in the limit 7 \ 0, the 
shear stress is obtained as a function of the temperature near T = Tg. The qualitative 
form of this function is displayed in the left side of figure |21 It is seen that the dynamic 
yield stress increases continuously from zero as T decreases from Tg. By contrast, when 
T = Tg and 7 > 0, we obtain 

CTxy ^ (16) 

with 5 = 3. This 5 corresponds to that appearing in the relation M ^ H^^^ describing 
critical phenomena for magnetic materials, and the value 5 = 3 is the mean field value 
for the Ising model. Note that equation (fT^ yields rj ^ The exponent 1 — 1/5 

is called the shear thinning exponent, and we find it to be 2/3 in this analysis. We 
also find that when T > Tg, a Newtonian regime appears near 7 = 0. This regime is 
connected to the power law regime that exists for higher shear rates. As seen in the 
inset of the left side of figure |2l this behavior corresponds to shear thinning. 

Next, we consider the case 63 > and find that here, a qualitatively different form of 
the stress is obtained. (We expect that there do exist models in which 63 > 0, although 
we believe that the model studied presently is not one of these.) One physically plausible 
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Figure 2. Shear stress in steady states in equation Here, the value plotted is 

that obtained in the limit 7 \ 0. The solid curve represents stable states, and the 
dashed line unstable states. The inset compares the dependence of axy on 7 in the 
three cases T>Ts,T = T^ and T < T^. Left : The case that 63 < 0. Right : The case 
that 63 > and 65 < 



situation is that in which 65 < 0, where we have the form displayed in the right side of 
figure 121 In this the temperature is decreased, the dynamic yield stress increases 

discontinuously from zero to a finite value at some temperature Tc(> Tg), as shown in 
the right side of figure |2l Furthermore, it is easily confirmed that shear thickening exists 
when T > Tc. 

5. Conclusion 

The main finding of this Letter is that the order parameter equation given in equation 
(I12|l provides a simple description of the nonlinear rheological phenomena exhibited by 
the systems considered here. We note that the order parameter in our model is defined 
to be the amplitude of the critical mode of the pair distribution function, which is 
related to the equal time correlation function of the density. This contrasts with the 
situation for the critical phenomena of liquid-gas phase transitions, for which the order 
parameter is defined in terms of the critical mode of the density. Consideration of this 
difference might help elucidate the essence of the glass transition. 

We conjecture that the coefficient 63 is negative for the model we study. This 
indicates that the dynamic yield stress appears continuously from at the transition 
point. This result is different from that obtained by the mode coupling theory that 
predicts the discontinuous appearance of the dynamic yield stress. Note that the recent 
result in Ref. ^1] might support the discontinuous onset, while the power-law behavior 
consistent with our result is still observed in the range 10"'^ < 7 < 10^^ in this report. 
Our simple theory in the present version might miss some important physical effects. 
We will study further in order to understand this discrepancy. 

Finally, we present remarks regarding the two assumptions we employed in our 
analysis. First, the perturbations to geq{r) are assumed to be restricted to the form 
h{r,t) = ip{r,t)lmY2^2{0,4>) for simplicity. Here, we note that we already performed 
the analysis without this restriction and obtained the same stability diagram [22] • In 
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this analysis, the five critical modes, which correspond to a^y, (Tyz-, cFzx-, <^xx — CTyy, and 
Czz ~ {<^xx + '^j/(/)/2, appear simultaneously at T = Tg as the result of the rotational 
symmetry. Then, five order parameters are defined in association with these critical 
modes. We will report the result in another paper. 

Second, with regard to the Kirkwood superposition approximation, we consider 
that it has some analogy with a mean field theory for critical phenomena. Therefore, 
we expect that the time evolution of the five order parameter fields (that depend on 
the spatial coordinate) under the infiuence of noise can describe rheological phenomena 
more precisely. This extension might modify the critical exponent 5, and also it might 
play an important role in the description of phenomena for the case T < Tg. We wish 
to develop such a theory as a natural extension of our analysis in this Letter. 
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